!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!>\brief
!!
!! Two stream instability
!!
!! \n
!!
!! For references, consider the following papers: 
!! <a
!! href="http://dx.doi.org/10.1016/S0010-4655(99)00247-7">[Nakamura,
!! Yabe, 1999]</a>, <a
!! href="http://dx.doi.org/10.1016/j.jcp.2009.11.007">[Crouseilles,
!! Mehrenberger, Sonnendr&uuml;cker, 2010]</a>.
!<----------------------------------------------------------------------
module mod_two_stream_inst_test

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_fu_manager, only: &
   mod_fu_manager_initialized, &
   new_file_unit

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_two_stream_inst_test_constructor, &
   mod_two_stream_inst_test_destructor,  &
   mod_two_stream_inst_test_initialized, &
   test_name,        &
   test_description, &
   coeff_init

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members
 character(len=*), parameter   :: test_name = "two_stream_inst"
 character(len=100), protected :: test_description(5)

 logical, protected :: mod_two_stream_inst_test_initialized = .false.

 real(wp), allocatable :: &
   k_wave(:),    & !< spatial wave number
   amplitude(:), & !< amplitude of the spatial perturbation
   vth(:),       & !< thermal velocity, main Gaussian
   vstream(:)      !< stream position in velocity space

 ! Input file ----------
 character(len=*), parameter :: &
   test_input_file_name = 'two_stream_inst.in'

 character(len=*), parameter :: &
   this_mod_name = 'mod_two_stream_inst_test'

 ! Input namelist, to be read from bump_on_tail.in
 namelist /input/ &
   k_wave, amplitude, vth, vstream

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_two_stream_inst_test_constructor(d)
  integer, intent(in) :: d(2)

  integer :: fu, ierr
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
          (mod_kinds_initialized.eqv..false.) .or. &
     (mod_fu_manager_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_two_stream_inst_test_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   ! Allocate the module variables
   allocate(k_wave(d(1)) , amplitude(d(1)) , vth(d(2)) , vstream(d(2)))

   ! Read input file
   call new_file_unit(fu,ierr)
   open(fu,file=trim(test_input_file_name), &
      status='old',action='read',form='formatted',iostat=ierr)
    if(ierr.ne.0) call error(this_sub_name,this_mod_name, &
      'Problems opening the input file')
    read(fu,input)
   close(fu,iostat=ierr)

   write(test_description(1),'(a,i1,a,i1)') &
     'Two-stream instability, space dimension ',d(1), &
                        ', velocity dimension ',d(2)
   write(test_description(2),'(a,*(e12.6,:," , "))') &  
     '       space wave numbers: ',k_wave
   write(test_description(3),'(a,*(e12.6,:," , "))') &  
     '  perturbation amplitudes: ',amplitude
   write(test_description(4),'(a,*(e12.6,:," , "))') &  
     '       thermal velocities: ',vth
   write(test_description(5),'(a,*(e12.6,:," , "))') &  
     '          stream position: ',vstream

   mod_two_stream_inst_test_initialized = .true.
 end subroutine mod_two_stream_inst_test_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_two_stream_inst_test_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_two_stream_inst_test_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   deallocate(k_wave , amplitude , vth , vstream)

   mod_two_stream_inst_test_initialized = .false.
 end subroutine mod_two_stream_inst_test_destructor

!-----------------------------------------------------------------------

 !> Two-Gaussian profile
 !!
 !! See \c mod_bump_on_tail_test for a very similar problem.
 pure function coeff_init(x,v) result(f0)
  real(wp), intent(in) :: x(:,:), v(:,:)
  real(wp) :: f0(size(x,2))

  real(wp), parameter :: pi = 3.1415926535897932384626433832795028_wp
  integer :: i
  real(wp) :: ncoeff, v2(size(v,2)), dx(size(x,2))

   ! first Gaussian distribution
   ncoeff = 1.0_wp/product(sqrt(2.0_wp*pi)*vth)
   v2 = 0.0_wp
   do i=1,size(v,1)
     v2 = v2 + ((v(i,:)-vstream(i))/vth(i))**2
   enddo
   f0 = 0.5_wp * ncoeff * exp( -0.5_wp*v2 )

   ! second Gaussian distribution
   ! ncoeff is the same
   v2 = 0.0_wp
   do i=1,size(v,1)
     v2 = v2 + ((v(i,:)+vstream(i))/vth(i))**2
   enddo
   f0 = f0 + 0.5_wp * ncoeff * exp( -0.5_wp*v2 )

   ! spatial perturbation
   dx = 1.0_wp
   do i=1,size(x,1)
     dx = dx + amplitude(i)*cos(k_wave(i)*x(i,:))
   enddo
   f0 = dx * f0
  
 end function coeff_init

!-----------------------------------------------------------------------

end module mod_two_stream_inst_test

